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^j- ■ Abstract 

o ■ 

' First- and second-order temperature driven transitions are studied, 

in a lattice gas driven by an oscillatory field. The short time dynamics 
study provides upper and lower bounds for the first-order transition 
points obtained using standard simulations. The difference between 
upper and lower bounds is a measure for the strength of the first- 
order transition and becomes negligible small for densities close to 
q one half. In addition, we give strong evidence on the existence of 

JD ■ multicritical points and a critical temperature gap, the latter induced 

by the anisotropy introduced by the driving field. 
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Far from equilibrium systems (FFES) are ubiquitous in nature and their 
theoretical understanding will contribute to the progress of scientific areas 
in physics, chemistry, biology, ecology, economy, etc. Since the theoretical 
development of non-equilibrium statistical mechanics is still in its infancy, 
a useful approach to FFES is to study simple models by means of various 
techniques such as numerical simulations, mean-field approximations, phe- 
nomenological scaling, field-theoretical developments, etc. Within the broad 
context of FFES, driven diffusive systems (DDS) Jl| have very recently re- 
ceived growing attention [§, [| |J; for reviews see e.g. || ||, The classical 
model for DDS was proposed by Katz et al. (KLS) [|1]] and is based on the 
equilibrium Ising model ||. Using the lattice gas language, the KLS model 
introduces an external driving field to the Ising model. However, due to this 
modification, the system now evolves towards a non-equilibrium stationary- 
state (NESS). In spite of considerable effort devoted to the study of the KLS 
model, there are still controversies on the understanding of numerical data 
||, |3], [| and its theoretical description is the subject of an ongoing debate 
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In this work we study a DDS subjected to the action of an oscillating 
driving field. One of the motivations for this approach is that the periodical 
field can be realized in numerous practical applications such as charged col- 
loids between the plates of a capacitor [jnj, electrophoresis experiments in 
pulsed fields |L3|], gas condensation in the presence of ultrasonic waves fl4 



segregation of granular materials in vibrating containers, etc. 

The aim of this work is to perform an extensive simulation study of the 
dependence of the temperature-driven transitions of the model on both the 
density of particles and the magnitude of the field. Measurements of station- 
ary properties combined to an study of the short time dynamics allow us to 
drawn a detailed phase diagram of the model that lead us to the discovery 
of a multicritical point. Furthermore, we developed a coupled mean field 
approach that yields results in agreement with the simulations. 

The model is defined on the square lattice assuming a rectangular geom- 
etry L x ,L y , using "brick wall" (periodic) boundary conditions across (along) 
the y— (x— )axis where the oscillatory field is applied, respectively. A lattice 
configuration r\ is specified by the set of occupation numbers n^j = {0, 1}, 
corresponding to each site of coordinates i.e. r\ = {riij}. Nearest- 

neighbor attraction with a coupling constant J > 0, is considered. So, in the 
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absence of a field the Hamiltonian H is given by 

H = -4J n i,3 n i',3'i (!) 

<ij,i'j'> 

where the summation is over nearest-neighbor sites only. The driving oscil- 
latory field E acts along the ±y— direction with half-period r. The coupling 
to a thermal bath at temperature T and the action of the field are considered 
through Metropolis jump rates, namely mm[l, exp — ({AH — eE(r)} /ksT), 
where ks is the Boltzmann constant, A7i is the change in 7i after the ex- 
change, and e = (—1,0,1) for a particle attempting to hop (against, or- 
thogonal, along) the driving field, respectively. For E = and half-filled 
lattices, the model reduces to the Ising model in absence of magnetic field. 
In the thermodynamic limit the Ising model exhibits a second-order phase 
transition at a temperature T/ = 2.2692.. J/ks- 

Monte Carlo simulations are performed on lattices of aspect ratios L x j L y = 
2 and 1, with 30 < L y < 480. T is reported in units of J/ks and E is given 
in units of J. The starting configuration is obtained by randomly filling the 
sample with probability p Q , which is also the density of particles that remains 
constant. One Monte Carlo time step (mcs) involves L x Ly trials. Data are 
obtained disregarding 10 6 mcs in order to allow the system to reach a NESS, 
and averages are taken over the subsequent 10 6 mcs. Using this procedure a 
single data point, as e.g. shown in figure 2, requires ~ 1 day of CPU time in 
an AMD 700 MHz processor. 

The model has also been studied by means of a coupled mean-field (CMF) 
approach. In order to write down the CMF equations the local density of 
particles p it j at site is defined which is the probability of finding a 

particle in this site. Due to normalization, one has p it j + h it j = 1, where hi j 
is the probability for the site to be empty. Then, one has to consider 
all events that may cause pij to change, pij may increase by the arrival 
of particles due to unbiased (biased) diffusion perpendicular (parallel) to the 
driving field, respectively. Similarly, the density may decrease by an outgoing 
flux of particles to neighboring sites. The implementation of the CMF leads 
to a set of L x Ly coupled non-linear differential equations. Here, we will only 
sketch out the form of such equations for the sake of space. Let r)'[(i, j); (?', f)\ 
be the configuration obtained from i] by interchanging the content of site 
with that of a neighboring site Then, the Metropolis rates are 

functions F of H(rf) - H(rj) - eE(r) = AH[(i,j); (i',f)] - eE(r). So, Pij 
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Figure 1: 3d plot of the density distribution characteristic of a NESS multi- 
stripped configuration obtained with the CMF method. L x = 80, L y = 40, 
T = 2.0, p = 0.50, E — 10 and r = 10. 

evolves in time according to: 

= hi,j{Pi+i,jF {AH[(i,j); (i + 1, j)],T} + pi- hj F{AH[(i,j); (i - l,j)),T} + (2) 

p i>j+1 F{AH[(i,j); + 1)],T,E(t)} + pij-xF{AH[(i, j); - l)],T,E(r)}} - 
Pid {h i+ld F{AH[(i,j); {i + l,j)],T} + h*- ld F{AH[(i,j)\ (i ~ hj)],T} + 
hi, j+1 F{AH[{i,j); + l)],T,E{r)} + h itj ^iF{AH[(i,j); - 1)], T, E{r)}}. 

The set of equations (3) is solved numerically starting from a random initial 
distribution of particles and using an integration time step of At = 0.25, in 
arbitrary units. Numerical integrations are performed until t = 25000 and 
averages are taken for t > 20000. In the CMF approach the excluded vol- 
ume interaction is taken into account in a probabilistic way and stochastic 
fluctuations are disregarded, in contrast to the Monte Carlo method which 
has intrinsic fluctuations and excluded volume is deterministically satisfied. 
However, the CMF approach is derived directly form the microscopies, so it 
contains the same symmetries than the lattice gas model. One advantage of 
the CMF method is that one can obtain the spatial mass distribution. In 
fact, figure 1 corresponds to a NESS where a multi stripped pattern is ob- 
served. An intriguing feature of driven dissipative systems is the occurrence 
of highly ordered and complex patterns as shown in figure 1. Since the system 
constantly gains (loss) energy from (to) the external field (thermal bath), re- 
spectively, the observed stationary states are by no means equilibrium states. 
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In fact, they are truly non-equilibrium steady states. 

In order to perform a quantitative investigation, the longitudinal order 
parameter (OP x ) is defined as the excess density, namely 

OP x = (RL x )- 1 ^ i \P(i)-p l (3) 

i=i 

where P(i) = {L y )~ l Y.% 

l riij is the density profile along the x— direction 
and R = (2p Q (l — p Q )) is a normalization constant. Similarly, OP y can also 
be defined. 

The dependence of the nature of the ordered phase on the period of the 
applied field has been investigated [T5[]. For temperatures below criticality, 
it is found that for short periods (say r < 4L y ) the system exhibits NESS 
with stripped patterns such as that shown in figure 1. However, for larger 
periods (say r > 4L y ) the system alternates between almost equilibrium 
states (AES) such as those corresponding to molecules in a gravitational 
field. The crossover from NESS to AES has a characteristic time of the 
order of r « 4L y . In this work, we are interested in the critical behavior 
of NESS so we have restricted ourselves to the case r = 10 mcs, without 
loosing generality because the same behavior will be valid for periods such 
as t < AL y for finite lattices and all periods in the thermodynamic limit. So, 
r plays an important role in this model. In fact, for the case treated in this 
work, namely r < 4L y , OP x is a well defined quantity independent of time t. 
However, for r > 4:L y , OP x and OP y are functions of time t, since the system 
alternates between AES as mentioned above. So, the half-period changes the 



nature of the problem and a crossover from NESS to AES is observed JTo 
In addition, since the oscillatory field causes the current of the driven gas 
averaged over long times to vanish, the symmetries of the model are different 
from those of the KLS model. From the theoretical point of view, this fact is 
essential to establish the universality class of the model, as will be discussed 
below. 

Figure 2(a) shows results obtained for E = 1. For low densities (p Q < 
0.15) the observed transitions are abrupt and exhibit strong metastability, so 
they are first-order. In contrast, for p Q > 0.40 one observes second-order or 
very weak first-order like behavior. Notice that for p Q = 0.20 and p Q = 0.40 
we have also included data which demonstrate the particle-hole exchange- 
invariance of the results. The existence of both first- and second-order tran- 
sitions can also be observed by using the CMF approach. These results are 
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Figure 2: (a) Plots of OP x versus T obtained for L x = 240, L y = 120, E = 1, 
r = 10 mcs and different values of p Q : o,p Q = 0.05; □, p = 0.075; A,p c = 
0.10; v,Po = 0.15; +,p D = 0.20; B,p = 0.80; A,p c = 0.40; ;p a = 0.60 and 
p = 0.50. The inset shows results obtained solving the CMF equations 
for E = 10 and r = 1. • p = 0.15, B,p = 0.30 o p G = 0.50. (b) Plots of 
OP x versus T obtained by using the Monte Carlo method for p Q = 0.05 and 
the values of the field indicated in the figure, (c) As in (b) but solving the 
CMF equations for p a = 0.30. 



in excellent agreement with Monte Carlo data, as shown in the inset of figure 
2(a). Figure 2(b) and 2(c) show that for low densities, T c depends on the 
amplitude of the field, so that the higher the field the lower the T c . Further- 
more, these figures also reveal that weaker first-order transitions are obtained 
for smaller amplitudes of the field. Remarkably, results obtained by means of 
the CMF approach exhibit the same decreasing trend than the Monte Carlo 
data. 

Figure 3 shows the phase diagrams obtained for a fixed lattice size {L X J L y = 
2, L y = 120) and two values of the driving field, namely E — 1 and E = 50 ~ 
oo. Using a method recently proposed for the study of the short time dy- 
namics of weak first-order transitions |16| it is possible to determine both 



lower and upper bounds for T c (p G ) valid in the thermodynamic limit and fur- 
ther generalize the phase diagram for E > 1. The idea behind the proposed 
method is based on the existence of two pseudo critical points at T* and T** 
near the weak first-order transition point T c with T* < T c < T**. These 
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Figure 3: Phase diagram, T c vs p Q , obtained from the data of figure 2(a). 
Empty (filled) symbols correspond to E = 1 {E = oo), respectively. On the 
left side, the symbols A and vy show the upper and lower bounds for T c as 
obtained by means of the short time dynamics study for E = oo. □ shows 
the location of the multicritical point. The full and the dashed curves, drawn 
on the right side, correspond to the best fit of the data obtained using eq. 
(|^). *, ■ and o show the location of the critical temperature of the Ising 
model T/ , the critical temperature predicted by eq. (H) for E — > 0, and the 
lower bound obtained for E = 0.01 using the short time dynamic analysis, 
respectively. 
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Figure 4: Log-log plots of the OP versus t, as obtained by means of the 
short time dynamics study. Averages are taking over 10 3 different runs using 
lattices of size L x /L y = 2,L y = 120. The following cases are shown: a) 
Starting from an ordered state, b) Starting from a fully disordered state, c) 
As in b) but for different lattice sizes, d) as in b) but for different driving 
fields. 



points can be obtained accurately from two short time dynamical processes 
starting from fully disordered and zero temperature states, respectively. In 
second-order transitions T* and T** overlap with the transition point T c , so 
the difference between T* and T** also gives a criterion for the weakness 



of the first-order transition [16]. Consider a system at T < T c (p ) and the 



evolution process from a fully disordered state. Due to the geometrical con- 
strained L x /Liy ^> 1 (0 « 0.2) |H] configurations at short times exhibit 
multi-stripped patterns that are long lived, only relaxing to the single stripe 
state after a time of the order t ~ L z x L y |TJ]]. Even in the case of square 
geometry, both the present model and the KLS model display multi-stripped 
configurations up to t ~ 10 5 mcs ||17||. It is then clear that the short time 



dynamics must be studied using an order parameter which takes into account 
multi-stripped configurations as that given by eq. (3). 

Our results for the short time dynamical behavior have been summarized 
in figure 4. For the used density (p = 0.16) power laws have been obtained 
for T** = 2.76 (figure 4(a)) and T* = 2.40 (figure 4(b)) starting from ordered 
and fully disordered states, respectively. Also, figure 4(c) shows that the 
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lower bound given by the short time dynamics is independent of the lattice 
size. Notice that the curves obtained for different aspect ratios are shifted 
but the power law behavior is obtained at the same temperature. The same 
results have been obtained for the upper bound, pointing out that the bounds 
drawn in the phase diagram (figure 3) are independent of the lattice size and 
consequently also valid in the thermodynamic limit. The transition points 
estimated using a finite lattice (figure 2) satisfy T* < T c < T** as would also 
do the true transition points in the L x ,L y — > oo limit. Also, the difference 
AT = T** — T* depends on the strength of the first-order transition while 
AT = at the second-order transition point for p a — 1/2. 

Coming back to the phase diagram, it is found that for p Q > 0.30, T C (E) 
steadily increases with the strength of the field, reaching a saturation value 
at T C (E = oo) ~ 1.41 T C (E = 0) for p Q = 1/2, in excellent agreement with 
results for the KLS model P, E|. However, for lower densities (e.g. for p Q < 
0.1, in figure 3) T C (E) steadily decreases when increasing the magnitude of 
the field. So, T C (E) exhibits opposite trends depending on the density and 
consequently, it is expected that for some characteristic density p^f (0.20 > 
Po 1 > 0.15) the critical temperature will be the same for all magnitudes of 
the driving field. Therefore, the point (p^ 1 ,T C (E, p^ 1 )) is a multicritical 
point, in the sense that for these special values of density and temperature 
this point is a critical point for all values of the amplitude of the field. Due 
to the observed symmetry, (1 — p^f ,T C (E, p^f)) is also a multicritical point. 

Assuming that the critical curves have the simplest form allowed by the 
symmetry of the system, we propose the following expression for the critical 
temperature: 

T c (E, Po ) =T c {^M2)-k 00 f{E){\±p'jyl? -k 00 {l-f{E)){\±p o y/P 1 £>0(4) 

where for E —>■ oo, A;^ is the coefficient of the higher order term and 
f(E) — » 0, respectively. Equation (||) can be thought as the first approx- 
imation to the phase coexistence curve, valid close to p a — 1/2, so that (3 is 
the order parameter critical exponent of the second-order transition. In order 
to fit eq. (U) to the data we will first summarize the symmetries present in 
our model. The model exhibits full translational and reflexion invariance as 
the Ising model, but the rotational symmetry is broken due to the anisotropy 
introduced by the field. If we consider short time scales, the up-down sym- 
metry is also broken by the field. However, a renormalization group study 
will consider the system at a coarse-grained level. Then, we expect that the 



9 



up-down symmetry will be restored at long time scales. Consequently, the 
present model displays the same symmetries than the randomly driven lattice 
gas with (3 = 1 jnj. Taking this value for (3, the critical curve for E = oo 
can be fitted using a single parameter, yielding koo = 15 ± 3. Assuming 
that f(E) = exp (-E), p* f is the only parameter left to be fitted, yielding 
p^ 1 = 0.160 ± 0.005 for E — 1 (see figure 4). Discrepancies between the fit 
and the data for densities far from p Q = 1/2 are expected since the expansion 
given by eq. (f|) holds close to that point only. Notice that equation (||) 
satisfies that (pf = 0.160 ± 0.005, T c {pf) = 2.59 ± 0.01)) is a multicritical 
point. This value is in agreement with the estimation performed using the 
short time dynamics study that gives T** = 2.76 > T c {p h J) > T* = 2.40 
(see figure 4). The existence of the multicritical point can also be confirmed 
by means of a short time dynamics simulations. In fact, figure 4(d) shows 
that plots of OP x versus t obtained for different fields (1 < E < oo) yield 
the same lower bound for T c {p^) given by T* = 2.40, independently of the 
strength of the field. This behavior is characteristic of the multicritical point, 
as observed in figure 3. It should be noticed that fits of the phase diagram 
assuming (3 = 1/2, as theoretically expected for the KLS model [0, are far 
from being satisfactory. Also, an excellent fit of the curve can be obtained 
assuming (3 = 1/4 (yielding (pf = 0.168 ±0.005, T c (pf) = 2.57±0.01)), but 
this value of the exponent is not supported by the symmetry considerations 
above mentioned. 

For the sake of comparison we have included in the phase diagram the 
critical temperature of the Ising model T* as well as the prediction of eq. 
dU) in the E — > limit. The latter is in excellent agreement with the lower 
bound estimate given by the short time dynamics method for E = 0.01. 
Notice that these estimations for the driven system are consistent with the 
location of the multicritical point that should also hold for E — > 0. These 
results show that, for p = 1/2, there is a gap in the critical temperature 
between the case E = (Ising model) and the limit E — > of the present 
model. Such a gap is expected to be even greater for p ^ 1/2 because in this 
case the coexistence temperature of the Ising model is lower than T* while 
the coexistence temperature of the driven diffusive system has a lower bound 
given by the multicritical temperature. The existence of these temperature 
gaps dramatically reflects the anisotropy introduced by the driving field and 
the non-equilibrium nature of the studied model. A physical explanation of 
this observation remains as an open question. 

In summary, the phase diagram of a DDS in the presence of an oscillatory 



10 



driving field is determined for E — 1 and E — oo. We give strong evidence on 
the existence of a multicritical point and a critical temperature gap separating 
the cases E = from E — > 0. To our best knowledge, these features have 
never been reported in the field of DDS. 
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